The present invention relates to a special purpose method and apparatus for efficiently computing forces or potentials in molecular dynamics.
In molecular dynamics, behavior of a liquid, a solid or a high polymer is considered to result from movements of constituent atoms and is investigated by simulating those movements with a computer. A certain proper region is taken from an actual substance, and a model involving coordinates, masses, charges, velocities, etc. of all the atoms included in that region is constructed on a computer. For example, in the case of a system (model) including 10,000 particles, a force acting on an ith particle is a sum of forces from the other 9,999 particles. Therefore, it is necessary to perform 9,999 force calculations and add up the calculated results. Since this operation is performed for every particle, 99,990,000 calculations are required in total. Based on the force thus determined, an acceleration of each particle is calculated and then its speed and coordinates after a lapse of a very short time are calculated. If determination of a new speed and coordinates of each particle is called "one-step progress of calculation," several thousand to several tens of thousand steps of calculation need to be performed in molecular dynamics.
As is apparent from the above, in general, in a system including N particles, the amount of calculation of forces is approximately proportional to N.sup.2 and the amount of calculation of speeds and coordinates is proportional to N. Therefore, the amount of calculation of forces becomes enormous for a large N, and it is difficult even for a supercomputer to conduct such calculation.
This is the reason why there are conventionally required special purpose computers for high-speed calculation of forces disclosed in: "FASTRUN: A Special Purpose, Hardwired Computer for Molecular Simulation," PTOTEINS: Structure, Function, and Genetics 11, pp. 242-253 (1991) and "GRAPE-2A: A Special Purpose Computer for Simulations of Many-Body Systems with Arbitrary Central Force," Proceedings of the 24th Hawaii International Conference on System Sciences, Vol. 1, pp. 171-180 (1992).
There are several kinds of forces which act between atoms. A van der Waals force acts between any types of atoms. If atoms have charges, a Coulomb force occurs between those atoms. If atoms are bonded to each other, there exists a bonding force that is in accordance with the bonding. Since in general the bonding force calculation is complex, it should be performed separately from the calculation of the other kinds of forces. Since only several atoms are connected to a single atom, the bonding force calculation is usually performed by a general purpose computer, not by a special purpose computer. That is, van der Waals forces and Coulomb forces are mainly calculated by a special purpose computer.
In many cases, the van der Waals force is expressed in the form of a well known L-J (Lennard-Jones) potential. This type of force F.sub.LJ is expressed as EQU F.sub.LJ =A(B.sup.6 r.sup.-7 -2B.sup.12 r.sup.-13) (1)
where r is a distance between particles concerned and A and B are constants determined by a kind of particles. Since the powers of r are negative large values, F.sub.LJ steeply decreases as r increases. Therefore, a cutoff distance r.sub.c is determined in accordance with the accuracy required for the calculation of F.sub.LJ. No calculations need to be performed for particle pairs whose r is larger than r.sub.c. This is illustrated in FIG. 2, in which the inside of a spherical surface indicated by a solid line is a system to be considered. A spherical surface having a radius r.sub.c is assumed around an ith particle indicated as a solid circle, and forces from only mesh-applied particles located inside the spherical surface need to be calculated (forces from particles indicated by void circles need not be calculated). Since the number of particles located around a single particle is considered to be substantially constant independently of the number N of particles included in the system, the amount of calculation of forces is proportional to N.
In view of the above, in the above-mentioned two conventional techniques, a particle pair list is prepared on a host computer and forces are calculated with a special purpose computer using the list. However, to prepare such a list, it is necessary to calculate all the distances between the ith particle and the other particles. Further, this operation needs to be performed for every particle, the total amount of calculation is proportional to N.sup.2 if the system includes N particles. In addition, since respective particles move, the particle pair list should be updated at proper time intervals. As a result, the proportion of the time for calculating the particle pair list to the total calculation time becomes not negligible as N increases.
Another subject to be considered when using the particle pair list relates to the capacity of a memory for storing the list. In the above-mentioned conventional technique (FASTRUN), the number of particle pairs is 583,267 when the number of particles included in the entire system is 2,204 and r.sub.c is 1.5 nm. If 32 bits (4 bytes) are required for the storage of one particle pair information, in this case about 2.3 megabytes are needed in total. Although the memory capacity of such an extent can be realized easily, the number of particles is on the order of 100,000 to deal with, for instance, a membrane of a living body. Further, if r.sub.c is increased to 2.8 nm to improve the calculation accuracy, in which case about 7,500 particles around a single particle should be considered, the capacity of about 1.5 gigabytes is required. It is very difficult to realize such a large capacity.
As described above, while the use of the particle pair list is convenient to calculate forces because the amount of calculation of forces is proportional to N, there exist problems that the amount of calculation to prepare the pair list is proportional to N.sup.2 and that the capacity of a memory for storing the list increases in proportion to N.
Another problem arises in dealing with a Coulomb force F.sub.C, which is expressed as EQU F.sub.C =CDr.sup.-2 ( 2)
where C and D are charges of particles concerned. In this case, since the power of r is small, F.sub.C does not decreases steeply as r increases. Since the Coulomb force is a strong, long-range force, it is not proper in essence to introduce the cutoff distance r.sub.c. In the above-mentioned two conventional techniques, the cutoff distance is used at the cost of the calculation accuracy.
On the other hand, since the van der Waals force and the Coulomb force are expressed by simple formulae (see Eqs. (1) and (2)), a special purpose computer can be constructed relatively easily by implementing calculation algorithms by hardware.
However, the bonding force cannot be calculated in a simple manner unlike the above case, but complex calculations need to be performed in accordance with the type and state of bonding. In a high polymer or the like, not only atoms directly bonded to a certain single atom but also atoms bonded to each directly bonded atom and farther atoms bonded to each atom that is bonded to a directly bonded atom need to be considered in the form of corrections. Since such corrections hardly processed by a special purpose computer and the number of atoms involved is several to ten odd, it is desirable to process those corrections on a general purpose computer (host computer) that controls a special purpose computer. Therefore, a certain means for separating particles to be subjected to calculations by the special purpose computer and particles to be subjected to calculations by the general purpose computer, as described below with reference to FIG. 3.
FIG. 3 schematically shows a state in which a certain substance is dissolved in water. In FIG. 3, a particle i indicated as a solid circle is an oxygen atom and mesh-applied particles j and k are hydrogen particles. These particles are bonded to each other to form a water molecule. An atom (or molecule) dissolved in water is denoted by 1. A spherical surface having the atom i at its center and having a proper short radius r.sub.a can be assumed so that only the atoms j and k are located inside the spherical surface (the other atoms are located outside it). Also in the case of a high polymer, it is possible to assume, by properly selecting r.sub.a, such a spherical surface as contains only the particles to be given special treatment. Further, it is in many cases proper to calculate by a host computer forces from particles that are not bonded to but very close to a particle of concern. Thus, it is very effective in the force calculation to separately calculate forces from particles inside and outside the spherical surface having the short radius r.sub.a.
The above-mentioned two conventional techniques have a function of calculating forces acting between a particular particle and all the remaining particles and automatically generating, at the same time, a list of numbers of particles in close proximity to the particular particle. However, the conventional apparatuses calculate forces from the particles of close proximity according to Eqs. (1) and (2) and includes those forces in the sum of forces. As described below, this will cause a problem in the calculation accuracy.
In the conventional technique (GRAPE-2A), a force is calculated by 32-bit floating point arithmetic (single precision). In the single precision type calculation, in which the mantissa portion consists of only 23 bits, the calculation accuracy of at most 7 orders (a little more than 6 orders) can be assured in decimal notation. Since the particles in close proximity to the particular particle impart very strong forces to the latter, if such forces are included in the sum of forces a portion that does not fit in the 23-bit mantissa portion is rounded. Therefore, even if corrections are later effected by the host computer based on the proximal particle list, highly accurate calculation results are not expected, because very large numbers (which should not be calculated according to Eq. (1)) having large errors have already been included in the sum.
Such a deterioration of the calculation accuracy could be prevented even in the conventional technique (GRAPE-2A) by properly modifying it. For instance, this could be done by controlling, with the host computer, the calculations of the forces from the particles in close proximity to the particular particle so that the constant A becomes 0 (constants A and B are determined by the kind of particles). To this end, however, the host computer needs to obtain in advance the proximal particle list by a certain measure and control, based on the list, the calculations so that A becomes 0. Since the proximal particle list should be changed for a different particular particle, the above operation needs to be performed every time the particular particle is switched.